{
 "metadata": {
  "name": "mass_spring_system-Copy0"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Mass Spring System:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this example, we are going to study a simple mass-spring system. This example serves to demonstrate the visualization handling abilities of mechanics-visualization module which is still in development process"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Dependencies:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For this system, following python modules are required:\n",
      "\n",
      "1) sympy=0.7.2\n",
      "\n",
      "2) numpy=1.6.2"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext autoreload\n",
      "%autoreload 2\n",
      "from sympy.physics.mechanics import *\n",
      "from sympy import symbols,cos,Eq,sqrt,sin\n",
      "import viz                       #This module is included with the notebook\n",
      "from IPython.display import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 49
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That is all we are going to need to carry out this visualization.\n",
      "\n",
      "Now, we will create symbols for coordinates, speed, force etc of the mass(bob)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q = dynamicsymbols('q')                  # Generalized coordinate\n",
      "u = dynamicsymbols('u')                  # Generalized speed\n",
      "f = dynamicsymbols('f')                  # Net Force applied to the box\n",
      "m = symbols('m')                         # Mass of body\n",
      "g, k, t, L = symbols('g k t L')               # Gravity,spring constant, time and length of unstretched spring\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 50
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we create an inertial reference frame, assign an origin, and set origin velocity to zero."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "I = ReferenceFrame('I')            # Inertial reference frame\n",
      "O = Point('O')                     # Origin point\n",
      "O.set_vel(I, 0)                    # Origin's velocity is zero"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we define the bob as a Particle .."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "P = Point('P')                      # Center Of Mass of the box\n",
      "P.set_pos(O, q * I.y)               # Set the position of P\n",
      "P.set_vel(I, u * I.y )              # Set the velocity of P\n",
      "bob = Particle('bob',P,m)   \n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 52
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We assign forces acting on the Particle.\n",
      "Coordinate axes are defined to be as : x-axis: running from left to right, y-axis: running from top to bottom, z-axis: running from behind the plane to above the plane.\n",
      "\n",
      "Following the above notion, Gravity is taken to be in (+I.y) direction, and force due to spring(F=Kx), in -I.y direction .."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = m*g* I.y - k*(q)*I.y    #Net force on bob.\n",
      "force = [(P,f)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 53
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Setting up equation of motion for the body:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "diff_eq = [q.diff(t) - u]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 54
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we call the Kane's Method to generate the equation of motion of the bob"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "kane = KanesMethod(I, q_ind=[q], u_ind=[u], kd_eqs=diff_eq) # Initialize the object\n",
      "fr, frstar = kane.kanes_equations(force, [bob])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 55
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The Differential equation of the motion of bob is given by the expression: fr + frstar = 0, i.e. fr = frstar"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fr, frstar"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 56,
       "text": [
        "([-g*m + k*q(t)], [m*Derivative(u(t), t)])"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, We need to solve this differential equation, for this we would be utilizing code generator sub-module in sympy mechanics, which is under development, and is not merged in the package.\n",
      "\n",
      "First packing the parameters in a list."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "parameters = [m, k, g, L]     #Mass, stiffness, gravity\n",
      "parameter_vals = [1, 1, 9.8, 1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 57
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we import the code generator module."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import code_gen\n",
      "right_hand_side = code_gen.numeric_right_hand_side(kane, parameters)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 58
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The numeric_right_hand_side method takes the Kane's object and parameters in symbolic form, and returns the differential equation in a numeric form, which can be integrated using odeint() method in SciPy.\n",
      "\n",
      "Now we integrate right_hand_side using odeint"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.integrate import odeint\n",
      "x0 = [1,0]        #Initial conditions, q=1,u=0.\n",
      "t = [i*0.1 for i in range(0,10)]    #Taking 10 time intervals of 0.1 sec\n",
      "numeric_vals = odeint(right_hand_side, x0, t, args=(parameter_vals,)) \n",
      "numeric_vals"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 65,
       "text": [
        "array([[ 1.        ,  0.        ],\n",
        "       [ 1.04396334,  0.87853405],\n",
        "       [ 1.17541411,  1.7482901 ],\n",
        "       [ 1.39303891,  2.6005778 ],\n",
        "       [ 1.69466326,  3.42688139],\n",
        "       [ 2.07727345,  4.21894471],\n",
        "       [ 2.53704657,  4.96885373],\n",
        "       [ 3.06938873,  5.66911561],\n",
        "       [ 3.66898093,  6.31273356],\n",
        "       [ 4.32983223,  6.89327675]])"
       ]
      }
     ],
     "prompt_number": 65
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The variable contains the positions and velocities, we need only positions for this 1D simulation, so we separate them, and pack them alongwith time."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Q = []\n",
      "for entries in numeric_vals:\n",
      "    Q.append(entries[0])\n",
      " \n",
      "simulation_coords = [Q,t]     #Packing positions and time in a single array"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 63
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Simulation:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For simulation and visualization, we are going to utilize, the Javascript rendering capabilities of IPython notebooks to render WebGL's.\n",
      "We are using Three.js Javascript library for 3D visualizations."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viz.visualize(simulation_coords)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "javascript": [
        "\n",
        "//Modifying AnimFrame for smart animations for different browsers\n",
        "window.requestAnimFrame = (function(){\n",
        "      return  window.requestAnimationFrame       || \n",
        "              window.webkitRequestAnimationFrame || \n",
        "              window.mozRequestAnimationFrame    || \n",
        "              window.oRequestAnimationFrame      || \n",
        "              window.msRequestAnimationFrame     || \n",
        "              function(/* function */ callback, /* DOMElement */ element){\n",
        "                window.setTimeout(callback, 1000 / 60);\n",
        "              };\n",
        "    })();\t\n",
        "    \n",
        "    \n",
        "    \n",
        "\n",
        "\n",
        "\n",
        "\n",
        "// fetching coorinates from Python output\n",
        "var coordinates = [[1.0, 1.043963344984979, 1.1754141133673666, 1.393038914946831, 1.6946632604568606, 2.0772734459187014, 2.5370465725650542, 3.0693887293725495, 3.6689809279475774, 4.3298322287782103, 5.0453396329483411, 5.8083540269903864, 6.6112516230633744, 7.4460101341520106, 8.3042889290908075, 9.1775123710715363, 10.05695550508608, 10.933831221694406, 11.799378083071952, 12.644947870758436, 13.462091881024476, 14.242645466951734, 14.978809604297208, 15.663228782059802, 16.289064520686569, 16.850063669767206, 17.340620914519064, 17.755834770309377, 18.091556555358487, 18.344431847834439, 18.511934002381899, 18.592389393806581, 18.584994138551849, 18.4898221240379, 18.307824279078766, 18.040819067522001, 17.69147431722892, 17.263280562295861, 16.760516169565314, 16.188204598152172, 15.552064198144958, 14.858451075166629, 14.114295582854174, 13.327033076139031, 12.504529621363243, 11.655003400373227, 10.786942600151775, 9.9090206348995338, 9.0300093260149481, 8.1586915020008028, 7.303773076883961, 6.4737961345647053, 5.6770535301135059, 4.9215060430446975, 4.2147028534381468, 3.56370610475483, 2.9750203422783539, 2.4545275208507893, 2.0074282314625376, 1.6381897415622193, 1.3505013604415244, 1.1472375770664169, 1.0304293357170493, 1.0012437446905, 1.059972416445093, 1.2060285540442111, 1.4379528151063226, 1.7534278887878063, 2.1493016519961436, 2.621618662853384, 3.1656596880458734, 3.775988848905417, 4.4465079390525402, 5.1705173525768693, 5.9407830266204043, 6.7496087214287481, 7.5889129180462751, 8.4503095585569792, 9.3251918617874043, 10.204818292734096, 11.080399920848826, 11.943188218052438, 12.784562497122385, 13.596116018183057, 14.369740008750897, 15.097704675556272, 15.772736431569612, 16.388090590239063, 16.937618724835971, 17.415830140490293, 17.817946708118495, 18.139950597083427, 18.378624470443459, 18.531583563738522, 18.597299557760884, 18.575115858989967, 18.465254090430427, 18.268811973271564, 17.987752296222297, 17.624883268536969], [0.0, 0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6000000000000001, 0.7000000000000001, 0.8, 0.9, 1.0, 1.1, 1.2000000000000002, 1.3, 1.4000000000000001, 1.5, 1.6, 1.7000000000000002, 1.8, 1.9000000000000001, 2.0, 2.1, 2.2, 2.3000000000000003, 2.4000000000000004, 2.5, 2.6, 2.7, 2.8000000000000003, 2.9000000000000004, 3.0, 3.1, 3.2, 3.3000000000000003, 3.4000000000000004, 3.5, 3.6, 3.7, 3.8000000000000003, 3.9000000000000004, 4.0, 4.1000000000000005, 4.2, 4.3, 4.4, 4.5, 4.6000000000000005, 4.7, 4.800000000000001, 4.9, 5.0, 5.1000000000000005, 5.2, 5.300000000000001, 5.4, 5.5, 5.6000000000000005, 5.7, 5.800000000000001, 5.9, 6.0, 6.1000000000000005, 6.2, 6.300000000000001, 6.4, 6.5, 6.6000000000000005, 6.7, 6.800000000000001, 6.9, 7.0, 7.1000000000000005, 7.2, 7.300000000000001, 7.4, 7.5, 7.6000000000000005, 7.7, 7.800000000000001, 7.9, 8.0, 8.1, 8.200000000000001, 8.3, 8.4, 8.5, 8.6, 8.700000000000001, 8.8, 8.9, 9.0, 9.1, 9.200000000000001, 9.3, 9.4, 9.5, 9.600000000000001, 9.700000000000001, 9.8, 9.9]];\n",
        "\n",
        "var link_length = coordinates[0][0];\n",
        "var time = coordinates[1][0];\n",
        "var i=0;\n",
        "//Setting other varaibles\n",
        "var controls,scene,camera,renderer,reset; \n",
        "var $reset_button;\n",
        "var $canvas\n",
        "var sphere_geometry;\n",
        "\n",
        "var sphere, link;\n",
        "var animation_request_id;\n",
        "\n",
        "\n",
        "//initial Position, dimensions of bob\n",
        "var radius = 0.1, segments = 100, rings = 100;\n",
        "\n",
        "init_canvas();\n",
        "init();\n",
        "visualize();\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\t\n",
        "function init_canvas(){\n",
        "\t\n",
        "\t//This function initiates a basic canvas, with trackball controls, \n",
        "\t//all drawing work occurs in init() function.\n",
        "\t\n",
        "\t// first of all , a renderer ...\n",
        "\trenderer = new THREE.WebGLRenderer();\n",
        "\t\n",
        "\t//show the IPython handle container\n",
        "\tcontainer.show();\n",
        "\t\n",
        "\t// create a canvas div\n",
        "\t$canvas = $(\"<div/>\").attr(\"id\", \"#canvas\");\n",
        "\t//giving background color\n",
        "\t\n",
        "\t$canvas.attr(\"style\",\"background-color:rgb(104,104,104)\");\n",
        "\t\n",
        "\t//For this particular, giving a top left, position time index\n",
        "\t\n",
        "\t$canvas.append(\"<p>click left and move mouse to rotate camera, hit reset button to reset camera</p>\");\n",
        "\t$canvas.append('<div id=\\\"pos_index\\\" style=\\\"margin-left:10px;margin-top:20px;\\\">T= ' + time + '<br /> Pos..='+(-link_length)+'</div>');\n",
        "\t\n",
        "\t// Adding our canvas to IPython UI\n",
        "\t\n",
        "\t\n",
        "\t// Now lets add a scene ..\n",
        "\t\n",
        "\t// set the scene size\n",
        "\t var WIDTH = 400,\n",
        "\t    HEIGHT = 300;\n",
        "\t    \n",
        "\tscene = new THREE.Scene();\n",
        "\t\n",
        "\t\n",
        "\t//Add a camera to the scene..\n",
        "\t\n",
        "\t// set some camera attributes\n",
        "\tvar VIEW_ANGLE = 45,\n",
        "\t    ASPECT = WIDTH / HEIGHT,\n",
        "\t    NEAR = 0.1,\n",
        "\t    FAR = 1000;\n",
        "\t    \n",
        "\tcamera = new THREE.PerspectiveCamera(  VIEW_ANGLE,\n",
        "\t                                ASPECT,\n",
        "\t                                NEAR,\n",
        "\t                                FAR  );\n",
        "\t    \n",
        "\t        \n",
        "\t// the camera starts at 0,0,0 so pull it back\n",
        "\tcamera.position.z = 10;\n",
        "\t\n",
        "\t\n",
        "\t// Add trackball controls\n",
        "\t\n",
        "\tcontrols = new THREE.TrackballControls(camera, renderer.domElement);\n",
        "\t\n",
        "\t\n",
        "\treset = function(){ controls.reset();}\n",
        "\t\n",
        "\t\n",
        "\tscene.add(camera);\n",
        "\t\n",
        "\t\n",
        "\n",
        "    $reset_button = $('<button/>').attr('style','margin-left:40px;').click(reset);\n",
        "    $reset_button.append('Reset Camera');                             \n",
        "    $canvas.append($reset_button)\n",
        "\t\n",
        "\t$canvas.append(renderer.domElement);\n",
        "\telement.append($canvas);\n",
        "\t// start the renderer\n",
        "\trenderer.setSize(WIDTH, HEIGHT);\n",
        "\t\n",
        "\t\n",
        "\t// Add  axes ...\n",
        "\t\n",
        "\tvar axesMaterial = new THREE.MeshLambertMaterial(\n",
        "\t{\n",
        "\t    color: 0xFFFFFF\n",
        "\t    \n",
        "\t});\n",
        "\tvar x_axis = new THREE.Mesh(\n",
        "\t   new THREE.CubeGeometry(WIDTH, 0.03, 0.03),\n",
        "\t   axesMaterial);\n",
        "\t   \n",
        "\tscene.add(x_axis);\n",
        "\t\n",
        "\tvar y_axis = new THREE.Mesh(\n",
        "\t   new THREE.CubeGeometry(0.03, HEIGHT, 0.03),\n",
        "\t   axesMaterial);\n",
        "\t   \n",
        "\tscene.add(y_axis);\n",
        "\t\n",
        "\tvar z_axis = new THREE.Mesh(\n",
        "\t   new THREE.CubeGeometry(0.03, 0.03, 10),\n",
        "\t   axesMaterial);\n",
        "\t   \n",
        "\tscene.add(z_axis);\n",
        "\t\n",
        "\t// create a point light\n",
        "\tvar pointLight = new THREE.PointLight( 0xFFFFFF );\n",
        "\n",
        "\t// set its position\n",
        "\tpointLight.position.x = 100;\n",
        "\tpointLight.position.y = 100;\n",
        "\tpointLight.position.z = 100;\n",
        "\n",
        "\t// add to the scene\n",
        "\tscene.add(pointLight);\n",
        "\t\n",
        "\trenderer.render(scene, camera);    \n",
        "\t\n",
        "\t\n",
        "}\n",
        "\t\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "function init(){\n",
        "\t\n",
        "    // create the sphere's material\n",
        "\tvar sphereMaterial = new THREE.MeshLambertMaterial(\n",
        "\t{\n",
        "\t    color: 0xCC0000\n",
        "\t});\n",
        "\n",
        "\t// create a new mesh with sphere geometry\n",
        "\t\n",
        "\tsphere_geometry = new THREE.SphereGeometry(radius, segments, rings)\n",
        "\tsphere = new THREE.Mesh(sphere_geometry\n",
        "\t   ,\n",
        "\t   sphereMaterial);\n",
        "    sphere.position.y = -link_length;\n",
        "    sphere.position.needsUpdate = true;\n",
        "    sphere.geometry.dynamic = true;\n",
        "\t\n",
        "\t// add the sphere to the scene\n",
        "\tscene.add(sphere);\n",
        "    \t\n",
        "    // draw!\n",
        "\trenderer.render(scene, camera);    \n",
        "    \n",
        "    }\n",
        "\n",
        "\n",
        "\n",
        "function visualize() {\n",
        "\t\n",
        "\t\n",
        "\tcontrols.update();\n",
        "\t\n",
        "\t\n",
        "    \n",
        "    renderer.render(scene, camera);\n",
        "\t\n",
        "\trequestAnimationFrame(visualize);    \n",
        "\t\n",
        "   \n",
        "}\n",
        "\n",
        "var animate  =function() { \n",
        "\t\n",
        "\t\n",
        "\t\n",
        "\tlink_length = coordinates[0][i];\n",
        "\ttime = coordinates[1][i];\n",
        "\tconsole.log('link_length = '+link_length)\n",
        "\tsphere.position.y = -link_length;\n",
        "\t\n",
        "\t\n",
        "\tconsole.log(sphere.position.y);\n",
        "\tconsole.log('radii = '+sphere.geometry.radius)\n",
        "\t$(\"#pos_index\").html('T= ' + time + '<br /> Pos..='+(-link_length));\n",
        "\t\n",
        "\trenderer.render(scene, camera);\n",
        "\ti+=1;\n",
        "\tif(i == coordinates[0].length)\n",
        "\t\ti=0;\n",
        "\t\t\n",
        "\t\n",
        "\tanimation_request_id = requestAnimationFrame(animate);    \n",
        "\t}\t\n",
        "\t\n",
        "var stop_animate  =function() { \n",
        "\t\n",
        "\t\n",
        "\tcancelAnimationFrame(animation_request_id);\n",
        "\t \n",
        "\t}\t\t\n",
        "\t\n",
        "\t\n",
        "var $animate_button = $('<button/>').attr('style','margin-left:40px;').click(animate);\n",
        "$animate_button.append('Animate');                             \n",
        "$canvas.append($animate_button)\n",
        "\n",
        "var $stop_animate_button = $('<button/>').attr('style','margin-left:40px;').click(stop_animate);\n",
        "$stop_animate_button.append('Stop Animation');                             \n",
        "$canvas.append($stop_animate_button)\n",
        "\n"
       ],
       "output_type": "display_data",
       "text": [
        "<IPython.core.display.Javascript at 0x45dcd90>"
       ]
      }
     ],
     "prompt_number": 64
    }
   ],
   "metadata": {}
  }
 ]
}